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ABSTRACT 

In earlier work we identified a global, non-axisymmetric instability associated with the presence of an 
extreme in the radial profile of the key function C(r) = (Efi/ n 2 )S 2 / r in a thin, inviscid, nonmagnetized 
accretion disk. Here, £(r) is the surface mass density of the disk, fi(r) the angular rotation rate, S(r) the 
specific entropy, T the adiabatic index, and n(r) the radial epicyclic frequency. The dispersion relation 
of the instability was shown to be similar to that of Rossby waves in planetary atmospheres. In this 
paper, we present the detailed linear theory of this Rossby wave instability and show that it exists for a 
wider range of conditions, specifically, for the case where there is a "jump" over some range of r in S(r) 
or in the pressure P{r). We elucidate the physical mechanism of this instability and its dependence on 
various parameters, including the magnitude of the "bump" or "jump," the azimuthal mode number, 
and the sound speed in the disk. We find large parameter range where the disk is stable to axisymmetric 
perturbations, but unstable to the non-axisymmetric Rossby waves. We find that growth rates of the 
Rossby wave instability can be high, ~ 0.2£!k for relative small "jumps" or "bumps" . We discuss possible 
conditions which can lead to this instability and the consequences of the instability. 

Subject headings: Accretion Disks — Hydrodynamics — Instabilities — Waves 



1. INTRODUCTION 



Glatzcl 1988|; |Goldrcich, Goodman, fc Narayan 1986 ; Kate 
1987; |Narayan, Goldrcich, fc Goodman 1987 ; Kojima et 
al. 19891; p'awley 1991|; [Zurck fc Benz 1986| ; see review by 
Narayan & Goodman 1989). In these studies, the "disk" 



Tho oontral problem of aoorotion diok thoory 



ldur 



standing tho mechanism of angular momentum transport , 



The angular momentum must flow outward in order that 
matter accrete onto the central gravitating object. Ear- 
lier work suggested hydrodynamic turbulence as the mech- 
anism of enhanced turbulent viscosity v t — ac s h, with 
a < 1 a dimensionless constant, c s the sound speed, and h 
the half-thickness of the disk (Shakura & Sunyaev 1973). 
However, the physical origin and level of the t urbulence 
is not established (see Papaloizou & Lin 1995 for a re- 
view). Rec e ntly, the magneto-ro tational insta bility ( Ve- 



likhov 1959 


Chandrasekhar 1960; 


Hawley 1998 


) has been studied in 



tions and shown to give a Maxwell stress sufficient to 
give significant outward transport of angular momentum 



( Balbus fc Hawley 1998| ); that is, a statistically averaged, 



effective Shakura-Sunyaev alpha parameter is found to be 
a ~ 0.01 (Brandenburg et al. 1995). However, significant 
questions remain. The simulation studies are local in that 
a shearing patch or box of the actual disk of size <C r is 
treated; and the boundary conditions on the top of the box 
have so far been unphysical. Furthermore, disks in some 
systems (e.g., protostellar systems) are predicted to have 
very small conductivity so that coupling between matter 
and magnetic field is negligible. 

Interesting fundamental questions still exist concerning 
pure hydrodynamic processes in accretion disks. One well- 
known hydrodyna mic instability of disks is the Papaloizou 



Pringlc instabil ity (Papaloizou & Pringle 1984; Papaloizou 



& Pringle 1985; hereafter PP) which has been extensively 



studie d both in l i near theory and by numerical simula- 



tions (Blaes 1985; Drury 1985; Frank fc Robertson 1988 



is taken to be a thin torus (finite height) or annulus (in- 
finite along z) with the radial width much less than the 
radius. The corotation radius, where the phase velocity of 
the wave u} r /m equals the angular velocity of the matter 
O(r), has a crucial role in that a wave propagating radially 
across it can be amplified. However, significant growth of a 
wave typically requires many passages through the corota- 
tion radius and this requires reflecting inner and/or outer 
boundaries of the disk. Thus the PP instability depends 
on inner and outer disk boundary conditions which are 
artificial for accretion disks. 

Recently, we pointed out a new R ossby wave instability 
of non-magnetized accretion disks ( Lovelace et al. 1999| : 
hereafter Paper I). Paper I shows that the local WKB 
dispersion relation for the unstable modes is closely anal- 
ogous to that for Rossby waves in planetary atmospheres. 
Rossby vortices associated with the waves are well known 
in planetary atmospheres and give rise for example to the 
Giant Red Spot on Jupiter (Sommcria, Meyers, & Swinney 
1988; Marcus 1989, 1990). 

In the cases considered in Paper I, the instability oc- 
curs when there is a "bump" in the radial variation of 
C(r) = (SO/k 2 )^ 2 / 1 ", where S(r) is the surface mass den- 
sity of the disk, f2(r) the angular rotation rate, S(r) the 
specific entropy, L the adiabatic index, and k(t) the ra- 
dial epicyclic frequency. Such a bump may arise from the 
nearly inviscid accretion of matter with finite specific an- 
gular momentum £ onto a compact star or black hole. The 
accreting matter tends to "pile up" at the centrifugal ra- 
dius r c — I 2 /(GM) (with M the mass of the central object) 
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Fig. 1. — Illustration of the nature of the Rossby wave instability for m = 3 for the case of a homentropic Gaussian 
bump centered at a radius r = 1. The gray-scale component of the figure represents the surface density perturbation <5£ 
with the darkest region the highest density. The contour lines and the arrows show the velocity field as seen in a reference 
frame rotating at the angular rate of the disk fl(ro). This velocity field is Av = <5v + <fi [rCl(r) — foO(ro)], where 5v is 
the usual Eulerian velocity perturbation. The figure shows that the perturbation of the disk involves three anticyclones 
of closed streamlines or "islands" which approximately coincide with the regions of enhanced densities <$£. The radial 
width of the "islands" has been exaggerated to make them visible. 



where it gives a radially localized bump in C(r). 

In accretion disks in some systems, such as active galac- 
tic nuclei, the weak self-gravity at large distances can be 
important in forming narrow rings (Shlosman, Begelman, 
& Frank 1990). Such rings would give a bump in C(r) of 
the form we have assumed. The bump in C(r) is crucial 
because it leads to trapping of the wave modes in a finite 
range of radii encompassing the corotation radius. Thus, 
reflecting inner and outer disk radii required for the PP 
instability are not required for the Rossby wave instabil- 
ity. In contrast with studies of the PP instability, we allow 
a general equation of state where the entropy of the disk 
matter can vary with radius. 

The present paper develops in greater detail the work 
of Paper I and presents the full linear theory analysis of 
the instability without some simplifying assumptions made 
there. The physical nature of the Rossby wave instability 
is shown in Figure |l| where we have overlaid the velocity 
field onto the surface density variations. The basic theory 
is given in §|| and §[| and results are given in In §[| 
we discuss the details of this instability and compare it 
with other hydrodynamic instabilities. Conclusions of this 
work are summarized in 80. 



2. EQUILIBRIUM DISKS 

We consider the stability of non-magnetized accretion 
disks. The disks are assumed geometrically thin with ver- 
tical half-thickness h <C r, where r is the radial distance. 
We use an inertial cylindrical (r, cj>, z) coordinate system. 
The surface mass density and vertically integrated pres- 



sure are: 



rh rh 

^( r ) = / dz p(r, z) , P{r) = I dz p(r, z) 



-h 



-h 



The equilibrium disk is stationary (d/dt — 0) and axisym- 
metric (d/d<f) = 0), with the flow velocity v w u<^(r)0; that 
is, the accretion velocity v r is assumed negligible. Self- 
gravity of the disk is assumed negligible so that $(r) = 
—GM/(r 2 + z 2 ) 1/2 , where M is the mass of the central 
object. But we will discuss the effects of self-gravity in 

For the axisymmetric equilibrium disk, the radial force 
balance is: 



r 



r n 2 



1 dP 



dr 



(1) 



The vertical hydrostatic equilibrium gives h (r) ~ {c s /v^)r, 
where 

P 



r- 



(2) 



the square of adiabatic sound speed and T is an effective 
adiabatic index discussed further in 

The focus of this paper is on the stability of equilibrium 
disks with a slowly varying background shear flow and a 
finite-amplitude density variation over a finite radial ex- 
tent. Specifically, we envision conditions where inflowing 
matter accumulates at some radius, for example, a cen- 
trifugal barrier. This accumulation may occur if matter 
is supplied at a rate exceeding the rate at which matter 
spreads by say a turbulent viscosity. Depending on the an- 
gular momentum distribution of the incoming flow, newly 
supplied matter can form a disk with a finite "bump" or 
"jump" in the surface density at say r with a radial width 



3 




Fig. 2. — Parameters of the initial equilibrium for the homentropic step jump case ('HSJ'): The pressure P(t)/Pq\ the 
corresponding angular velocity Q(r)/£l K ; the square of the epicyclic frequency K 2 (r)/r2|-; and the effective sound speed 
Cs/vk- The pressure jump is derived from a jump in surface density £(r) with an amplitude A = 0.65 over Ar/ro = 0.05. 
This results in a pressure jump (1 +^4) r w 2.3 for T = 5/3. The disk rotation is close to Keplerian except near the region 
of the pressure jump. Note that k 2 is everywhere positive. 



Ar. When the disk surface density is sufficiently small, 
the disk optical depth is small compared with unity and 
cooling by radiation is efficient. As more matter is ac- 
cumulated, the disk becomes optically thick so that heat 
builds up inside the disk and can be vertically confined, 
then pressure forces start to have an important role for the 
disk dynamics. The present work is directed at optically 
thick disks with significant pressure forces. 

We model the mentioned "bumps" and "jumps" with 
fairly general radial profiles of £(r) and P{r). These pro- 
files are used to obtain the corresponding Sl{r) profile using 
equation (Q). Specifically, we consider two surface den- 
sity distributions: One is has a "step jump" from to 
£2 > Si over Ar and this jump is surrounded by a 
smooth background flow, 



S 



1 



A 
2 



tanh 



ro 



At 



1 



(3) 



where = £o( r Ao) - ' 3 is the surface density for the back- 
ground disk, So is its value at ro and (3 characterizes the 
slope (it is —3/4 in a standard a disk model). Quantities 
A and Ar measure the amplitude and width of the jump 
respectively, and ro is radius of the jump. The second case 
we consider is a "Gaussian bump," 



— = 1 + (A - 1) exp 



ro 



Ar 



(4) 



we can further derive or specify the radial distributions 
of temperature T(r)/To and pressure P(t)/Pq based on 
£(r), where T and Pq are the values at r for the back- 
ground disk. We also set vo = r Cl to be the Keplerian 
speed at ro, and define the dimensionless sound speed as 
Cq = T(Pq/T,q)/vq. In most of the following analysis, we 
take r = I, Ar/r = 0.05, T = 5/3, and c /t>o = 0.1. 



2.1. Homentropic Disks 

For homentropic flows, we can determine P/Pq and 
T/T using E/S alone because of the relations P/Pq = 
(E/E ) r and T/T = (P/P )/(£/£ ). Based on these, we 
present two sample initial configurations: one is a step 
jump with A = 0.65 which we call case 'HSJ' (home- 
ntropic step jump); the other is a Gaussian bump with 
A = 1.25 which we call case 'HGB' (homentropic Gaus- 
sian bump). The adiabatic index is 7 = 5/3 and the 
width Ar/ro = 0.05. Figures || and || show the profiles 
of P/P(r ), Q(r)/n K (r), K 2 (r)/tf K (r), and c s /v K (r Q ) for 
the 'HSJ' and 'HGB' cases. 

Here, we have chosen the parameter (3 to be zero in the 
smooth profile of £*. We have studied the dependence of 
our results on (3 and find that for < (3 < 3/4, our results 
are essentially independent of (3. Thus, we omit further 
discussion of the background disk. 



2.2. N on- Homentropic Disks 

For nonhomentropic flows, we have to specify T/Tq in 
addition to S/So m order to determine P/Pq and conse- 
quently r^(r). Such conditions were studied in Paper I for 

4 In actual fluids, the pressure depends on say both density S and entropy S. In this paper, we use the term homentropic to indicate that 
the pressure depends only on density with the entropy a constant, P = P(S). The term barotropic is sometimes used for the same situation. 
Similarly, we use nonhomentropic (instead of nonbarotropic) to describe a flow where the entropy is not a constant and P = P(T,,S). 



where again, A and Ar measure the height and width of 
the bump, respectively. 

We assume an ideal gas equation of state, P oc ET. 
Depending on whether the flow is homentropic 4 or not, 
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Fig. 3. — Similar to Figure ^]but for the case of an homentropic Gaussian bump ('HGB'). The peak of the surface density 
is A = 1.25 at r — r , which gives a peak of the pressure A T ~ 1.45 for T — 5/3. Again, k 2 is everywhere positive. 
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Fig. 4. — Similar to Figure but for the case of a nonhomentropic step jump ('NSJ'). The jump amplitude of the surface 
density and temperature is A — 0.52, which gives a pressure jump of (1 + A) 2 ~ 2.3. Again, k 2 is everywhere positive. 



a Gaussian bump, which we will not repeat here. Instead, 
we study a case where there a simultaneous step jump in 
both surface density and temperature. In order to avoid 
too many parameters, we will assume that the width and 
amplitude of the jumps are the same for surface density 
and temperature, which are described by equation (^). We 
will call this initial configuration as case 'NSJ' (nonhome- 
ntropic step jump). Figure || shows the dependences of 
different variables with A = 0.52 and Ar/ro = 0.05. 

As we show later, the nonhomentropic cases have some 
quantitative differences from the homentropic cases but 
the essential physics is the same. 

3. PERTURBATIONS OF DISK 



We consider small perturbations to the inviscid Euler 
equations. The perturbations are considered to be in the 
plane of the disk. Thus the perturbed surface mass density 
is: E = E + <5E(r, (j), t); the perturbed vertically integrated 
pressure is: P = P + SP(r, <fi,t); and the perturbed flow 
velocity is: v = v + <5v(r, 0, t), with <5v = (Sv r , Sv^, 0). 
The equations for the 2D compressible disk are: 



~Dt 

m _ 
~Dt ~ 



E V • v = 



-Ivp. 

E 



(5) 



(0) 
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D 
Dt 



P 



= 0, 



(7) 



where 



D d 

— = h v • V 

Dt dt 



(see Paper I). Here, S = P/S r is referred to as the entropy 
of the disk matter. Equation ([?]) corresponds to the adia- 
batic motion of the disk matter. The rigorous relation for 
the perturbation of the disk is of course D(p/ p 7t ) / Dt = 0, 
where p, p and 7t are the 3D pressure, density and adi- 
abatic index, respectively. (We reserve the symbol 7 for 
mode growth rate.) For slow perturbations of the disk 
matter (time scales < h/c s ~ 1/^), we have h oc y/T so 
that T oc £2(7t-i)/(7t+i) j anc i consequently P oc S r with 
T = (37t-l)/(7t+l). Equations (|) - (0) are the vertically 
integrated 2D Eulcr equations. The equations assume that 
the disk thickness does not change rapidly with radius or 
time, \dh/dr\ <C 1 and \dh/dt\ <C Qr. The above equa- 
tions are closed with an equation of state. The ideal gas 
law is applicable for the conditions of interest P = ST///, 
where p is the mean mass per particle. Notice that the 
pressure is in general a function of both the surface den- 
sity S and the entropy S (or equivalently, the temperature 
T). This is general in contrast with the commonly made 
assumption that P = P(S). 

Following the steps of Paper I, we linearize the equa- 
tions by considering perturbations oc f(r) exp(/m0 — itut), 
where m = ±1,±2, etc. is the azimuthal mode number 
and u = iv r + ij is the mode frequency. We use ^ = SP/T, 
as our key variable which is analogous to enthalpy of the 
flow. (In a homentropic flow considered by PP, \P is the 
enthalpy.) The basic equations are given in Paper I, but 
for completeness we also give the main equations here: 



iAtoSY, = V ■ (Stfv) , 
which is from the continuity equation (|^); 



S . S<5w 
SY, = — $ + 1 



Y>5v r = i T 



Auj 



AujL 



— I -2A-..VV 



(8) 

(9) 
(10) 



S(5«0 = T 



-kr. 



Auj 

IT 



AujQL s L p 



\ T k2 







(-£)] 


) 2il 2 





which are from Euler's equation 



Here, 



(11) 



m 
r 



K 

Aiv{r) 
Hr) 

Ts 



1 d{r 4 n 2 ) 

r 3 fir ' 
= u> — mf2(r) = uj r — mf2(r) + ?7 , 
= S n/[ K 2 - Alu 2 - c 2 /{L s L p )} , 



V 

r 



AlnS 

dr 



d ( P 

dr 



(12) 



where |Ts|, \L S \ and \L V \ are the radial length scales of 
the surface density, entropy and pressure variations, re- 
spectively. They are related as 



1 

L,, 



1 

IT 



1 

Ts 



(13) 



Furthermore, we can obtain 



rT 

7T 



s* 

„2 



2k A r 



T 1 / rT V 4£u J" 





OL 2 . ' r \ VLL S ) ~ r Awi s ^ Auj 2 flL s L p 

(see Paper I), which can in turn be written as 
+ B(r)*' + C(r)# = , 

where 



* , (14) 



(15) 



and 



Cl = k\ 



B(r) 
C{r) 

k 2 - Auj 2 



1 T' Cl' 
r + ~T ~ Tl 

-ci - c 2 , 



ft T' 



(16) 
(17) 



(18) 



(>2 



l-L' s B{r)+Ak^/Auj k 2 c 2 JAoj 2 -l 



L 2 



L s Lp 



(19) 



For homentropic flow, L s — > 00. Thus the coefficients in 
the above equations simplify to give 



s n 



k* — Alu* 
C(r) = -ci . 



(20) 
(21) 



In this limit, our equation ( jig ) is the same as that given 
previously (for example, PP)7^ 

For any given equilibrium disk, we can answer the fol- 
lowing questions: (1) Are there unstable modes with pos- 
itive growth rates 7 > 0? (2) What is the nature of the 
radial wave functions \t (r)? (3) What is the dependence of 
to r and 7 on the initial equilibrium? (4) And, what is the 
physical mechanism(s) of the instability. Equation (15) 
allows the determination of 10 — uj r + ij and identification 
of mode structure for general disk flows, both stable or 
unstable. 

3.1. Axisymmetric Stability 

It is important to know if the equilibrium disk including 
the "jump" or "bump" is stable to axisymmetric pertur- 
bations. Rather than solving the non-local axisymmet- 
ric stability problem (which can be obtained from equa- 
tion 15), we consider the local stability criterion for ax- 
isymmetric perturbations. This condition is expected to 
be a sufficient condition for non-local stability. If the 
disk pressure is neglected, then the Rayleigh criterion 
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'NSJ' case 




Radius (r ) 



Fig. 5. — A comparison of k 2 (upper panel) and A 2 (lower panel) for the nonhomentropic step jump ( 'NSJ', see Figure 
||). Note the difference in magnitude of the two frequencies. The disk is locally stable according to the Solberg-Hoiland 
criterion, n 2 + N 2 > 0. 



k 2 = (1 /r 3 )d(il 2 r 4: ) / dr > implies axisymmetric stability 
(Drazin & Reid 1981, ch. 3; Binney & Tremaine 1987, ch. 
6). The more general condition including the pressure is 
the Solberg-Hoiland criterion, 

K 2 (r) + N 2 { r) > , where N 2 = -- - --) 

(22) 

(Endal & Sofia 1978). Here, A is the radial Brunt- Vaisala 
frequency for the disk due to the radial entropy variation. 
For homentropic flow, N 2 is zero. Typically N 2 (r) is neg- 
ative for standard disk models (Shakura & Sunyaev 1973). 
Although this is an unstable situation for convection in the 
absence of rotation, it is stabilized by the rotation since 
| A'' 2 1 is usually much smaller than n 2 in an approximate 
Keplerian disk. Figure || shows illustrative profiles of k 2 (r) 
and N 2 (r) for the case 'NSJ' with A = 0.52. The quantity 
k 2 + N 2 remains positive. 

As the amplitude of the bump /jump A increases, we 
can find a critical value where condition given by equation 
(§§) is violated. This happens at Arit ~ 0.74,0.55,1.26 
for cases 'HSJ', 'NSJ' and 'HGB', respectively. Physically, 
this means that the local specific angular momentum pro- 
file is perturbed enough that it liberates energy when lo- 
cally exchanging two fluids radially. Actual axisymmetric 
instability is expected to occur only at values of A larger 
than these A c values. Note also that the critical values 
are dependent upon some other parameters, such as Ar, T 
and whether the flow is homentropic or not. The impor- 
tant point we note here is that we find non-axisymmetric 
instability at values of A appreciable less than the values 
Ac ■ 

3.2. Methods of Solving Equation ( |7j[ ) 

Since the eigenfrequency lu is in general complex, equa- 
tion (16) is a second-order differential equation with com- 
plex coefficients which are functions of r. If we discretize 



equation ([l5| ) on a grid i — 1, ... N and combine it with ap- 
propriate boundary conditions, the problem becomes that 
of finding the complex roots of the determinant of an Ax A 
tridiagonal matrix. Upon finding the roots, we can then 
solve for the corresponding complex wave function \I/(r). 

We use a Nyquist method to find the uj eigenvalues. 
For this we integrate along a closed path in the complex 
uj— plane and thereby find regions where roots (and poles) 
reside. We then use a Newton root-finder to locate the 
roots accurately. Once a root is found for a particular set 
of parameters, we can find other roots by varying the pa- 
rameters slowly. There are two apparent singularities in 
equation (|l5|), one is at the corotation resonance where 
Alu = 0, which occurs only on the real uj axis. The 
other is where W = n 2 — Alu 2 — c 2 /(L s L p ) — 0, which 
is a generalized form of the Lindblad resonance condi- 
tion for L s finite. Goldreich et al. (1986) showed that 
the first singularity is a real singularity, whereas the sec- 
ond is spurious and is actually a regular point. They 
considered the case when L s — > oo, but their conclusion 
still holds for L s finite. Thus, extra caution is needed in 
searching for roots in the complex plane with 7 > 0, be- 
cause the condition W = can be satisfied both when 
7 = and Alo 2 = k 2 — c 2 /(L s L p ), or when Au> r = and 
7 2 = — k 2 + c 2 /{L s L p ). The latter case is generally only 
possible when n 2 — c 2 /{L 8 L p ) < 0, which can be true if 
the pressure bump /jump is strong enough. 

3.3. Boundary Conditions 

The inner and outer boundary on ^(r) in equation (16) 
are important. In most previous studies where a torus is 
considered, the boundary condition is that the Lagrangian 
pressure perturbation vanishes at the free surfaces (both 
above and below the torus and at its radial limits) . In the 
present problem, the "bump" or "jump" is embedded in 
a background disk so that the appropriate boundaries are 
necessarily at a large distance from the bump but still in 
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the disk. Inspection of the function C in equation (15) 
shows that r/r < 0.8 and r/r > 1.2 correspond to large 
distances from the bump. In these outer regions we can 
use a WKB representation of with <3> = /(r)exp{is}. 
We then have 



is 



„'2 



iBs' + C = 



(23) 



where a prime denotes a derivative with respect to r. Away 
from the bump, the s" term is small and it can be dropped. 
Thus, we find for r <C Tq, 



1 



0' 

'77 



mf2 

c s 



15 1 

16 r 



(24) 



Thus, toward the inner boundary, 

r -15/16 



exp{i J s'dr}| oc 

increasing as r decreases. This is indeed the de- 
pendence found for the calculated wavefunctions discussed 
in later sections. Similar procedure can be applied to the 
outer boundary too. 

The physical meaning of this radiative boundary condi- 
tion can be understood by studying the dispersion relation 
obtained from the WKB approximation at the inner and 
outer parts of the disk. Let k 2 ~ s' 2 « TZe(C), we have 



(u> r — mQ) 2 = k 2 + k 2 c 2 



(25) 



Since k 2 w Q 2 at both n and r<i away from the perturbed 
region, the dispersion relation is then simply 



mQ ± y/Q 2 + k 2 c 2 s 



(26) 



As it turns out that the most unstable mode usually 
has w r close to corotation at ro, i.e., mQ(ro), the above 
equation implies that we need to choose "-" at r\ since 
uj r <C mfi(ri) and similarly "+" at outer boundary. Fur- 
thermore, the group velocity v g = dw/dk r w ^fc s k r /\k r \ 
at inner and outer radius, respectively. Since physically 
we expect the wave to propagate away from the central 
region, this requires that v g < at r\ and v g > at r2, 
which means k r > at both r\ and r%. 

The computational results presented here we take the 
inner boundary to be at r\ = 0.6ro and the outer bound- 
ary at ri — 1.4ro, where ro is the radius of the "bump" or 
"jump." We find that there is essentially no dependence 
of our results on the values of r\ and 

We have experimented with other types of boundary 
conditions such as vanishing pressure or vanishing radial 
velocity. The overall properties of the unstable modes 
show very little change, but the relative phase between 
the real and imaginary parts of \& of course changes. 

3.4. Initial Value Problem 

A different approach to solving the linearized equations 
(|| - is to treat them as an initial value problem. We let 



/ Ul 



( <5£/£o 
Sv r /v 

\ SP/Po 



(27) 



where all variables have the same meaning as in §0. As- 
suming that the <j) dependence of u has the form exp(zm^>), 



we can rewrite equations ([s] - as a system of first order 
PDE in one dimension, i.e., 



du 
~dt 



G(r) 



du 
Or 



H{r) :u = 



(28) 



where it is straightforward to write down the elements of 
matrices G and H. It can be shown that these equations 
are hyperbolic. G or H can be diagonalizcd individually 
but not simultaneously. 

We have written a code to solve equations (28) using 
a 2nd order MacCormack method. The idea is to deter- 
mine the growth rates of different initial perturbations. Of 
course, solutions of (28) will eventually be dominated by 
the most unstable mode. From this solution we can get 
the growth rate, the real frequency, and the radial wave 
function for all 4 variables. 

However, we recommend the method described in § [3.2| 
because it allows a comprehensive search for unstable 
modes and the corresponding eigenfunctions. The initial 
value code for equation (28) can then be used to verify 
the behavior of these modes. We omit a detailed compar- 
ison of the results obtained from both types of codes but 
merely report that we have found excellent agreement be- 
tween them on growth rates, mode frequencies, and radial 
structure of eigenmodes. 

4. RESULTS 

We present the solutions to equation ( [l5J)_ for different 
types of initial equilibrium as discussed inSH. In Paper I, 
we showed that a necessary condition for instability is that 
the key function L{r) = (Ef2/ n 2 )S 2 ^ r have an extreme as 
a function of r. This condition is a generalized Raylcigh 
inflexion point theorem for compressible and nonhomen- 
tropic flow. All three initial equilibria presented earlier 
satisfy this necessary condition for instability. Here we 
discuss in detail the behavior of these unstable modes. 

4.1. Step Jump Case 'HSJ' with A = 0.65 

The initial equilibrium for this case is shown in Figure 
£| For m = 3, we find that the most unstable mode for 
this equilibrium has a growth rate 7/fio ~ 0.154 and a 
real frequency w r /(mfio) ~ 0.92. 

It is informative to plot the effective "potential," C(r) 
of equation JIe|). We can neglect B(r) because its magni- 
tude is much smaller than C(r). Thus, equation ( |l5| ) can 
be simplified to give 



C(r)* = 



(29) 



Then the function C(r) is analogous to the quantity 
E — V{r) in quantum mechanics except that here C(r) is 
complex (Paper I). The real and imaginary parts of C(r) 
are shown in Figure ^| (the upper panel). First, we consider 
the region of 0.95 < r/r < 1.05, which is most affected by 
the presence of the surface density jump. There is a neg- 
ative real part of — C which is analogous to the "potential 
well" in quantum mechanics where bound states are pos- 
sible, though again the fact that C is complex prevents 
a direct analogy. The two sharp "peaks" in this region 
are not singularities but result from the two extrema of 
1/(k 2 — Aw 2 ). In the region r/r < 0.8, the function — C 
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Fig. 6. — The effective "potential well" — C(r) in equation (fi~3) for the most unstable mode with to = 3, found using 
parameters of cases 'HSJ', 'NSJ', and 'HGB'. The solid and dashed curves are the real and imaginary parts of the function 
— C(r), respectively. The unstable modes are excited at radii near ro inside the potential well. The unstable modes are 
trapped in this potential well. This trapping, however, is not absolute since there is a finite "probability" for modes to 
tunnel through into both the inner and outer parts of the disk. 



is dominated by c\ ~ —(m£l/c s ) 2 oc — r~ 3 , whereas for 
r/r > 1.2, — C is dominated by c\ w ~{uj r /c s ) 2 oc — 
const. Now consider an unstable mode that is excited 
in the "potential well" around ro. The positive poten- 
tial around r/ro = 0.9 and 1.1 causes this mode to be 



evanescent in this region. The potential is negative for 
r/ro < 0.8 and r/ro > 1.2 so that there will be a finite 
probability for this mode to "tunnel" through the poten- 
tial "barriers". In other words, a mode excited around ro 
by the surface density jump can "radiate away" into both 
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Fig. 7. — Radial eigenfunction for the perturbed surface density, for the radial and azimuthal velocity perturbations, and 
for the pressure perturbation, for the most unstable mode of the homentropic step jump ('HSJ') for m = 3. The solid, 
dashed and long-dashed curves are the amplitude, the real, and the imaginary parts of the eigenfunctions, respectively. 
All wave functions show their largest amplitudes near the deepest position of the potential well (slightly away from ro). 
The amplitudes decrease going away from ro. The wave-like oscillations towards the inner boundary is due to the radiative 
boundary condition. Similar oscillations also occur towards the outer boundary, but their amplitude is too small to be 
evident. These features of the wave functions are consistent with the potential well shown in the upper panel of Figure pi 
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Fig. 8. — Similar to Figure ^ except that parameters for the homentropic Gaussian bump ('HGB') are used. Again, the 
structure of the eigenfunctions is consistent with the potential well shown in the bottom panel of Figure ^|. For example, 
the two peaks in the pressure eigenfunction corresponds to the two minima of the potential. 



the inner and outer parts of the disk. This is indeed the 
case as seen in Figure [7] which shows the radial eigenfunc- 
tions of various physical quantities derived from 'J'(r) us- 
ing equations (|9|)-(|ll|). In obtaining these eigenfunctions, 
we have used outward propagating sound wave boundary 
(i.e., radiative) conditions discussed earlier. The relative 
phase shift between real and imaginary parts indicates this 
propagation. Note that the "wavelength" of the unstable 
modes is at least 2Ar > h; that is, the 2D approximation 
is satisfied. Furthermore, we have maintained the actual 
relative amplitude among all 4 variables. 

4.2. Step Jump Case 'NSJ' with A = 0.52 

For 77i == 3, the most unstable mode for this nonhomen- 
tropic equilibrium has a growth rate 7/^0 ~ 0.176, and a 
real frequency ui r / '(mf^o) ~ 0.93. These values are slightly 
different from the homentropic case presented above with 
a slightly higher growth rate. Note that the two cases have 
roughly the same jump in total pressure but the homen- 
tropic case requires a higher surface density jump. This 
will always be the case if T < 2. Consequently, 'NSJ' case 
has a slightly deeper "trap" which is shown in both the 
k 2 and the effective potential — C as depicted in Figure || 
(the middle panel). The eigenfunctions for all 4 variables 
are almost the same as in the case 'HSJ' (cf. Figure g), 
therefore we do not show them here. 

The comparison between cases 'NSJ' and 'HSJ' confirms 
one of the conclusions in Paper I that the variation of the 
temperature is more effective than that of the surface den- 
sity in driving this instability. 

4.3. Gaussian Bump Case 'HGB' with A = 1.25 

For 777, = 3, the most unstable mode for this case oc- 
curs with a growth rate 7/^0 ~ 0.21, and a real frequency 
oj r / (mflo) ~ 0.98. The effective potential — C for this case 
is shown in Figure [3] (the bottom panel) . The correspond- 
ing eigenfunctions are shown in Figure @. 



It is interesting to note that for the rather small Gaus- 
sian bump used here (A = 1-25), the instability has a 
higher growth rate than the step jump cases considered 
above. This difference can be understood as due to the 
deeper "trap" produced by the Gaussian bump. This can 
be seen by comparing the k 2 /Vl 2 K profiles for both cases 
using Figure || and |[ Since d(il/ilfc)/dr is positive at 
both the rising and declining edges of a Gaussian bump, 
and is positive only at the rising edge of a step jump, 
the k 2 profile has two regions that are higher than Vl 2 K 
in 'HGB' compared to just one such region in 'HSJ'. In 
general, larger k 2 indicates stronger stability, and they 
directly translate into the "forbidden" regions in the po- 
tential structures shown in Figure ^ (the bottom panel). 
Thus, a mode excited around r$ is well trapped by the 
"walls" at r/r » 0.8 - 0.9 and 1.1 - 1.2. The 'HSJ' case 
also has a positive spike at t/tq ~ 1.05, but it does not 
provide good trapping because the spike is too narrow. 
Again, we have used the outward propagating sound wave 
boundary conditions discussed earlier so that there is a fi- 
nite probability for an unstable mode to tunnel through 
the potential barriers. 



4.4. Instability Threshold and Maximum Growth Rates 

Consider now the dependence of the growth rate and 
mode frequency on the bump /jump amplitude A and the 
azimuthal mode number m. As the amplitude A decreases, 
the growth rate of the instability is expected to decrease. 
At small enough A > 1, the instability should turn off (Pa- 
per I). It is clearly of interest to determine the minimum 
A for instability. 

Figure O shows the growth rate and mode frequency as a 
function of A for the three cases considered, all for m = 5. 
The threshold values are .Athrc « 0.19,0.16, and 1.08, for 
cases 'HSJ', 'NSJ', and 'HGB', respectively. Note that the 
values of *4thre depend on Ar because the radial gradient 
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Fig. 9. — Dependences of the mode frequencies and growth rates on the amplitude of the surface density jump/bump A 
with in = 5 for the three cases considered. The vanishing of the growth rate for A < Athres indicates the thresholds for 
individual cases. The vertical dashed lines show the critical val ues of A where k 2 + N 2 = 0. For larger values of A the 
flow violates the condition for local axisymmetric stability (see §3.1). 




Azimuthal Mode Number m 



Fig. 10. — Dependences of the mode frequency uj r and growth rate 7 on the azimuthal mode number to. The filled dots, 
squares and triangles are for 'HSJ', 'NSJ' and 'HGB' cases, respectively. 



of the surface density S(r) is important in giving instabil- 
ity. Roughly speaking, a factor of ~ 1.2 "jump" or an 8% 
"bump" in the surface density are sufficient to cause the 
Rossby wave instability. 

The vertical dashed lines in Figure ^ indicate the A val- 
ues beyond which part of the flow has k 2 < as discussed 
In the step jump cases, the the Rossby wave in- 



3.1 



stability continues smoothly through the point where k 2 
changes sign. In case 'HGB', the instability also exists for 
values of A where k 2 < in part of the flow. There is 
a continuous increase in growth rate as A increases (not 
shown here), although there seems to be a "kink" at the 
point where kappa 2 starts to have a range of negative val- 



ues. We do not pursue here the situation when both the 
Rossby and the axisymmetric instabilities are present si- 
multaneously. We emphasize that the Rossby instability 
discussed here has a substantial growth rate (< O^f^K^o)) 
when the flow is stable to axisymmetric perturbations. 

Figure [n] shows the dependences of the growth rate and 
mode frequency on the azimuthal mode number m for the 
three cases. The peak of the growth rates around to = 4, 5 
probably results from a preferred azimuthal wavelength 
in comparison with the radial wavelength of the unstable 
modes. 
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Fig. 11. — Dependences of the mode frequency uo r and growth rate 7 on the sound speed for the homentropic step jump 
('HSJ') for m = 5. 



4.5. Dependence on Sound Speed 

Because the Rossby instability depends critically on the 
pressure forces, its growth rate has a strong dependence 
on the magnitude of sound speed. This is shown in Fig- 
ure [n] for the 'HSJ' case for m — 5. As the sound speed 
decreases, the pressure forces are no longer strong enough 
to perturb the rotational flow, the instability disappears. 
This is similar to the threshold seen in Figure || when the 
pressure jump/bump is too small. Judging from the plot, 
the growth rate increases roughly linearly with c s . 

Note that the threshold value of c s also depends on Ar. 
As c s decreases, so does the thickness of the disk. Thus, 
the allowed Ar can be smaller so long as the wavelength 
of the unstable modes is larger than the disk thickness. 
This is a necessary condition for 2D approximation to be 
satisfied. 



5.1. 



5. DISCUSSION 

Origin of Initial Equilibrium 

The condition for the Rossby wave instability (here- 
after, RWI) discussed here can be understood in terms 
of the Rayleigh's inflexion-point theorem (Drazin & Reid 
1981, p. 81), which gives a necessary condition for instabil- 
ity. If the accretion disk is close to Keplerian everywhere 
with temperature and density being smooth powerlaws, 
the RWI does not occur. The local extreme value of the 
key f unction C(i) considciud hcic is caused by Lhc local 



step j ump or Gauooian bump in ourfacc donoity (or poooi — 
bly in temperature if the flow is nonhomcntropic) . How- 
ever, the considered profiles of E(r) and P(r) are not the 
only ones which may lead to instability. For example, a 
profile with a local extreme in potential vorticity distri- 
bution may also give instability. An important aspect of 
the RWI is that the required threshold for instability is 
quite small. Typically, a 10 — 20% variation of S(r) over 
a lengthscale slightly larger than the thickness of the disk 
gives instability. 



The question is, then, can such an initial equilibrium 
exist in real astrophysical systems? The answer to this 
question requires detailed knowledge about how matter is 
initially "brought" in towards the gravitating object. One 
plausible situation is that accreting matter is gradually 
stored at some large radial distance. Subsequently, after 
say a sufficient build up of matter then efficient accre- 
tion can proceed. This is perhaps the physical situation 
in the case of low mass X-ray binaries where the systems 
go through episodes of outbursts with long, quiescent in- 
tervals between them (see Tanaka & Shibazaki 1996 for 
a review). In protostellar disks, it is possible that dif- 
ferent regions of the disk have different coupling strength 
between matter and the magnetic field. This could lead 
to accumulation of matter at large distances where the 
disk is non-magnetized. At smaller distances, the accretion 
may be due to the Maxwell stress from due to turbulent 
magnetic fields arising from the Balbus-Hawley instability 
(Brandenburg 1998). Matter accumulation over some fi- 
nite extent in radius could most likely lead to enhanced 
gradients of surface density and/or temperature, thus sat- 
isfying the conditions for the RWI discussed here. 

Another aspect of the problem is the role of vertical con- 
vection. We have mostly discussed our instability in the 
2D limit where the vertical variation of physical quantities 
is all averaged away. In real disks, however, there are cer- 
tain processes in the vertical direction that cou ld occur on 
a fast times c ale, such as vertical convection (Papaloizou 



k Lin 1995| ; |Klahr, Hcnning, fc Kley 1999| ). The verti- 



cal convection is perhaps not a main concern here since 
it helps to bring the matter in a vertical column to the 
same entropy, which is assumed in our present study. For 
the horizontal motion, as we emphasized before, the initial 
equilibria we have studied are all stable to the local convec- 
tive instabilities. So, the RWI will always be the dominant 
instability unless there are other instabilities (in the {r, 4>} 
plane) that have lower thresholds and faster growth rates. 
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5.2. Influence of Self- Gravity 

For accretion disks in some systems, for example, active 
galactic nuclei and protostellar systems, self-gravity may 
be important at large distances and lead to axisymmet- 
ric instability. The WKB dispersion relation for to = 
perturbations of a self-gravitating disk is 

lu 2 = n 2 - 27rGS|fc r | + k 2 c 2 (30) 

(Safronov 1960; Toomre 1964). The least stable, smallest 
uj 2 mode occurs for k r = irGY,/c 2 , and for this wavenum- 
ber, uj 2 — n 2 (\ — l/Q). Thus, the disk is unstable if 
Toomre's parameter Q = nc s / '(itGYj) < 1. The thresh- 
old for instability Q = 1 corresponds to a surface mass 
density in the case of an AGN disk of 



J thre 



irG 



70.7- 



c s /v K \ ( M 



cm 2 V 0.01 ) \ 10 8 M ( 



Ipc 



(31) 

where we assumed k w Qk (Shlosman, Begelman, & Frank 
1990). 

Near the threshold with Q somewhat less than unity, 
this instability will act to make axisymmetric rings with 
radial wavelength A r 2'kc s ISIk = 2iTh. These rings have 
a form similar to the case of our "Gaussian bump." Only a 
relatively small linear growth of the rings is needed for the 
non-axisymmetric RWI to set in with a very large growth 
rate. It is important to note that the RWI is capable of 
transporting angular momentum and thereby facilitating 
accretion whereas the axisymmetric instability is not. 

5.3. Minimum Surface Density and Relevance to X-ray 

Novae 

We mentioned earlier the important role of the pressure 
(oc c 2 ) for the RWI. This leads to an important physical re- 
quirement on the minimum surface mass density E c of the 
disk for the RWI to be important. Heat must be confined 
in the region of RWI for a time much longer than rota- 
tion period in order that the mode grows without radia- 
tive cooling. This confinement requires a minimum surface 
density which depends on the opacity and specific heat of 
matter. This translates to roughly S c ~ 10 2-3 g cm~ 2 . 
In other words, if S is too small, the RWI will not oc- 
cur and there will be no Rossby vortex induced accretion. 
But once enough matter has accumulated so that S ex- 
ceeds this threshold S c , RWI sets in and there is efficient 
accretion due to Rossby vortices. This may correspond to 
the distinctly different states of accretion in X-ray binary 
systems where sources are often observed to be "active" 
or "quiescent". 

As shown by Tanaka & Shibazaki (1996), X-ray no- 
vae (especially black hole candidates) show long quies- 
cent intervals in X-rays despite that the companion star 
is still continuously feeding mass to the compact object at 
~ 10 15 ~ 16 gm s -1 . This strongly implies a mass accumula- 
tion at some large distances without much accretion going 
on. The accumulated mass over an interval of ~ 50 yrs 
will be ~ 1.5 x io 24 ~ 25 gm, which gives a surface density 
of £ - 5 x 10 2 ~ 3 g cm~ 2 with a size of - 3 x 10 10 cm. This 
critical X is interestingly close to the value we discussed 
above. So, we believe it is worthwhile to pursue RWI as 
an alternative mechanism for causing an outburst in X-ray 
novae, in addition to the usual disk instability models. 



5.4. Comparison with the Related Work 

It is difficult to make a direct comparison of the RWI 
with the original Papaloizou & Pringle instability. In PP 
instability, a torus and/or an annulus has two edges and 
the existence of certain unstable modes critically depends 
on these edges (i.e., the so-called "principal branch"). In 
our study, the surface density bump case can perhaps be 
viewed as a torus except that it is embedded in a back- 
ground Keplerian shear flow. Here, we discuss two major 
differences between the two instabilities. One is the treat- 
ment of boundary conditions. In our study, the propagat- 
ing sound wave boundary condition provides a natural and 
important link of the unstable region with the surrounding 
flow. This allows us to directly apply such instability to 
thin Keplerian accretion disk of a much larger radial ex- 
tent. In the nonlinear regime, this "leakage" will allow the 
unstable modes to grow further nonlinearly and impact a 
large part of the disk flow also. The second difference is the 
fact that, unlike PP instability, the rotational profile f2(r) 
is not taken as a single power-law in our case. Instead, 
enhanced epicyclic frequency occurs in association with 
the surface density jump/bump boundaries. This natu- 
rally creates the "potential well" that allows the unstable 
modes to grow and provides the trapping at the bound- 
aries, as illustrated by the structure of the "potential" in 
Figure ^|. The corotation radius occurs within this poten- 
tial well. 



5.5. Nonhomentropic vs. Homentropic 

In Paper I we have already given a necessary condition 
for instability with respect to 2D non-axisymmetric dis- 
turbances, that is, for a key function 



C(r) = T{r)S 2/T {r) 



(32) 



to have a local extreme, where = z • (V x v)/S is 
the potential vorticity, S — P/S r is the entropy. This is 
a generalized Rayleigh criterion w hich was originally de - 
rived for incompressible flows (cf. prazin fc Rcid 1981 ). 
Most of studies on PP instability have assumed homen- 
tropic flow where entropy of the whole fluid system is a 
constant. The inclusion of S(r) in equation (32) makes it 
perhaps more applicable to real astrophysical disks where 
entropy is usually not a constant. As we have shown in this 
paper, the onset of the RWI, however, does not depend on 
having an entropy gradient. On the other hand, as shown 
in Paper I and by the comparison of 'HSJ' and 'NSJ' cases, 
temperature variation is more effective in causing the in- 
stability. Furthermore, entropy gradient introduces more 
features into the problem, most notably that the poten- 
tial vorticity of the flow is no longer conserved due to the 
net thermodynamic driving of vorticity from the fact that 
VT x VS 7^ (see Paper I). 



5.6. Implications for Angular Momentum Transport 

In order to evaluate the angular momentum transport 
from this instability, it is necessary to evaluate terms which 
include second order terms. However, the present linear 
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Fig. 12. — The Reynolds stress (arbitrary units) at radius r as derived from the velocity variations in linear theory. The 
positive values of this quantity indicates outward transport of angular momentum. 



theory is valid only to first order. Instead of the angu- 
lar momentum transport, we calculate the Reynolds stress 
derived from the velocity perturbations < SvrSv^ >. This 
term is expected to be important for the angular momen- 
tum transport. The Reynolds stress at each radius can be 
calculated as, by averaging over 0, 

f 27r dCl 
Jo dr 

= — nr 2 T, ^— \lZe(5v r )lZe(5vA) + Im(5v r )Im(5vA,)] . 
dr 

_ (33) 
The result is shown in Figure \L2 for the 'HGB' case with 
m = 5 (results for 'HSJ' and 'N5J' cases are similar). This 
instability indeed causes a mainly positive Reynolds stress 
which corresponds to an outward transport. 

A full evaluation of the angular momentum and matter 
transport due to the RWI evidently requires non-linear hy- 
drodynamic simulations. We have begun to perform these 
simulations and will present our findings in a forthcoming 
paper (Li, Lovelace, & Colgate 1999). Briefly, we have 
observed the production of largescale 2D vortices in the 
{r, 4>} plane and these vortices are observed to survive the 

shear flow for many revolutions of the disk (~10). There is 
also indication of outward angular momentum transport, 
which shows promises for this instability being a robust 
mechanism of angular momentum transport. 

One important point we want to emphasize is that the 
transport process regulated by these large scale vortices 
is inviscid, highly dynamic and nonlocal. This is funda- 
mentally different from the standard a disk models where 
the flow is assumed to be viscous, hence the transport is 
at small scale (i.e., local) and a quasi-stationary state can 
usually be reached. In fact, the notion of a statistically 
stationary a value and stationary accretion may be an in- 
correct physical picture for the accretion in some highly 
variable systems such as X-ray binaries. 



6. CONCLUSIONS 

We have developed a detailed linear theory of the Rossby 
wave instability associated with an axisymmetric, local 
surface density jump or bump in a thin accretion disk. 
The instability is termed Rossby due to its WKB disper- 
sion relation which is analogous to that for Rossby waves 
in planetary atmospheres. Rossby vortices associated with 
the waves are well known in planetary atmospheres and 
give rise for example to the Giant Red Spot on Jupiter 
(Sommeria, Meyers, & Swinney 1988; Marcus 1989, 1990). 
The flow is made unstable due to the existence of local ex- 
treme value of a generalized potential vorticity C(r) which 
includes the radial variation of entropy. Depending on the 
parameters, the unstable modes are found to have sub- 
stantial growth rates ~ 0.2f2(ro), where ro is the location 
of enhanced surface density gradient. These modes are ca- 
pable of transporting angular momentum outward. Since 
this instability relies on pressure forces perturbing the ro- 
tational flow, it requires a minimum sound speed, or a 
local "heat" content. We expect that the disk must be 
optically thick in order for the instability to be important. 

It is important to understand the nonlinear interactions 
between this instability and the rest of the disk, and to 
see whether the instability is effective in transporting an- 
gular momentum globally. An important aspect of the 
present instability is that the unstable modes propagate 
into the surrounding stable disk flow, thus allowing the 
instability to impact a large part of the disk. We have 
performed preliminary nonlinear hydro simulations of this 
instability and will report the detailed results in a forth- 
coming publication. The simulations have confirmed the 
the linear growth of the Rossby wave instability and they 
have shown the vortices have long lifetimes (many orbital 
periods). 
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Kato and C. Fryer. HL gratefully acknowledges the sup- 
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